
#include <navierstokes.hpp>
#include <ns_cl.cpp>

#include <functionSup.cpp>

#define STABILISATION 0

#define WEAKCL 0

namespace Feel
{
namespace FSI
{

    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    //NavierStokes<Dim, Order, Entity>::updatePtFixe(element_fluid_type & U,
    //element_fluid_type & Uold,
    //element_fluid_velocity_type & mesh_velocity,
    //e_map_element_type & depl )
    NavierStokes<Dim, Order, Entity>::updatePtFixe( element_fluid_type & Uold, sparse_matrix_ptrtype& A , vector_ptrtype& F)
    {
#if 0
        using namespace Feel::vf;

        std::cout << "PtFixe : Assemblage start\n";

        //space_fluid_ptrtype XhBis = U.functionSpace();

        //vector_ptrtype F( M_backend->newVector( Xh ) );
        //sparse_matrix_ptrtype A( M_backend->newMatrix(Xh,Xh) );

        auto mesh = this->mesh();
        auto Xh = this->functionSpace();


        element_fluid_type V( Xh, "v" );
        element_fluid_0_type u = U_fluid->template element<0>();
        element_fluid_0_type v = V.template element<0>();
        element_fluid_1_type p = U_fluid->template element<1>();
        element_fluid_1_type q = V.template element<1>();

        element_fluid_0_type uold = Uold.template element<0>();
        element_fluid_1_type pold = Uold.template element<1>();

        //--------------------------------------------------------------------------------------------------//

        double EE = 21*1e5;
        double sigma = 0.28;
        double mu = EE/(2*(1+sigma));
        double lambda = EE*sigma/((1+sigma)*(1-2*sigma));

        auto Iddd = oneX()*trans(oneX()) + oneY()*trans(oneY());
        auto Fv = Iddd + gradv(*M_aux);
        auto Ev = 0.5*(gradv(*M_aux)+trans(gradv(*M_aux)) ) + 0.5*trans(gradv(*M_aux))*gradv(*M_aux);
        auto Sv = lambda*trace(Ev)*Iddd + 2*mu*Ev;

        //--------------------------------------------------------------------------------------------------//

        //Tenseur des deformations
        auto deft = 0.5*(gradt(u)+trans(gradt(u)) );
        //        auto def = 0.5*(grad(v)+trans(grad(v)) );

        //Matrice Identitee
        auto Id = oneX()*trans(oneX()) + oneY()*trans(oneY());
        // Tenseur des contraintes (trial)
        auto Sigmat = -idt(p)*Id + 2*M_mu*deft;

        // Tenseur des contraintes (test)
        //        auto Sigma = -id(q)*Id + 2*mu*def;

        //force volumique :
        //auto f = -rho*cst(9.81)*vec(cst(0.),cst(1.));
        auto f = vec(cst(0.),cst(0.));
        // Vitesse d'entree du fluide
        auto Vit_Entree = 10*Vitesse_Entree*( Py()*(1.0-Py()) )*oneX();

        //--------------------------------------------------------------------------------------------------//

        form2( Xh,Xh, A, _init=true );

        bool M_isNLTerme=true;
        if (M_isNLTerme)
            {
                //terme non lineaire
                form2( Xh,Xh, A )  +=
                    integrate (elements(mesh), _Q<3*Order>(),
                               rho*trans(id(v))*(gradt(u))*(idv(uold) - idv(*M_mesh_velocity)) );
            }

        // sigma : grad(v) sur Omega
        form2( Xh, Xh, A ) +=
            integrate( elements(mesh), _Q<2*(Order)>(),
                       + trace(Sigmat*trans(grad(v))) );

        // terme assurant l'incompressibilite
        form2( Xh, Xh, A ) +=
            integrate( elements(mesh),// _Q<2*(Order-1)>(),
                       + divt(u)*id(q) );


#if 1
        // terme necessaire pour condition forte avec petsc
        form2( Xh, Xh, A ) +=
            integrate( elements( mesh ),// _Q<3*uOrder-1+3*(GeoOrder-1)>(),
                       +0*div(v)*idt(p)+0*divt(u)*id(q) +0*idt(p)*id(q) );
#endif
        /*
        //terme en plus pour div=0 (voir these goncalo)
        form2( Xh, Xh, A ) +=
        integrate( elements( mesh ),// _Q<3*uOrder-1+3*(GeoOrder-1)>(),
        divv(uold)*trans(idt(u))*id(v) );
        */

        // forces volumique
        form1( Xh, F, _init=true ) =
            integrate( elements(mesh), _Q<Order>(),
                       trans(f)*id(v) );


    if (false)
        {
            form2( Xh,Xh, A ) +=
                integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                           //- trans( (-idv(p)*Id + 2*mu*defv)*N())*id(v) );
                           - trans( (/*-vf::cst(10)*Id*/ + 2*M_mu*deft)*N())*id(v) );
        }

            form1( Xh, F ) +=
                integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                           //- trans( (-idv(p)*Id + 2*mu*defv)*N())*id(v) );
                           + trans( (-vf::cst(1e4)*Id)*N())*id(v) );


        //bool M_isMoveDomain=false;
        if (M_isMoveDomain)
            {
                //imposition de la contrainte
                form1( Xh, F ) +=
                    integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<3*Order>(),
                               + trans(Fv*Sv*N())*id(v) );
                form1( Xh, F ) +=
                    integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<3*Order>(),
                               + trans(Fv*Sv*N())*id(v) );
            }
        else if (M_weakCL)
            {
#if 0
                //imposition de la contrainte
                form2( Xh, Xh, A ) +=
                    integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<3*Order>(),
                               - trans(Sigmat*N())*id(v) );
                form2( Xh, Xh, A ) +=
                    integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<3*Order>(),
                               - trans(Sigmat*N())*id(v) );
#endif
            }
        //--------------------------------------------------------------------------------------------------//
#if WEAKCL
        if (M_weakCL)
            {
                // terme de bord classique
                form2( Xh, Xh, A ) +=
                    integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                               - trans(Sigmat*N())*id(v) );
                //+ trans(idv(uold)*trans(idt(u))*N())*id(v) );

#if 0
                form2( Xh, Xh, A ) +=
                    integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
                               - trans(Sigmat*N())*id(v) );
                //+ trans(idv(uold)*trans(idt(u))*N())*id(v) );
#endif
                // pas de terme car contrainte imposee
                /*
                  form2( Xh, Xh, A ) +=
                  integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                  - trans(Sigmat*N())*id(v) );
                  //+ trans(idv(uold)*trans(idt(u))*N())*id(v) );
                  form2( Xh, Xh, A ) +=
                  integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                  - trans(Sigmat*N())*id(v) );
                  //+ trans(idv(uold)*trans(idt(u))*N())*id(v) );
                  */

                //terme de bord weak boundary
                form2( Xh, Xh, A ) +=
                    integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                               - trans(Sigma*N())*idt(u) // symetrisation
                               //+ trans(idv(uold)*trans(id(v))*N())*idt(u)
                               + penalbc*trans(idt(u))*id(v)/hFace() );
#if 0
                form2( Xh, Xh, A ) +=
                    integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
                               - trans(Sigma*N())*idt(u) //symetrisation
                               //+ trans(idv(uold)*trans(id(v))*N())*idt(u)
                               + penalbc*trans(idt(u))*id(v)/hFace() );
#endif
                form2( Xh, Xh, A ) +=
                    integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                               //- trans(Sigma*N())*idt(u) //pas de symetrisation si contrainte imposee!!
                               //+ trans(idv(uold)*trans(id(v))*N())*idt(u)
                               + penalbc*trans(idt(u))*id(v)/hFace() );
                form2( Xh, Xh, A ) +=
                    integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                               //- trans(Sigma*N())*idt(u) //pas de symetrisation si contrainte imposee!!
                               //+ trans(idv(uold)*trans(id(v))*N())*idt(u)
                               + penalbc*trans(idt(u))*id(v)/hFace() );


                form1( Xh, F ) +=
                    integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                               - trans(Sigma*N())*Vit_Entree // symetrisation
                               //+ trans(idv(uold)*trans(id(v))*N())*Vit_Entree
                               + penalbc*trans(Vit_Entree)*id(v)/hFace() );

                form1( Xh, F ) +=
                    integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                               //  - trans(Sigma*N())*idv(mesh_velocity) //pas de symetrisation si contrainte imposee!!
                               //+ trans(idv(uold)*trans(id(v))*N())*idv(mesh_velocity)
                               + penalbc*trans(idv(mesh_velocity))*id(v)/hFace() );
                form1( Xh, F ) +=
                    integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                               //- trans(Sigma*N())*idv(mesh_velocity) //pas de symetrisation si contrainte imposee!!
                               //+ trans(idv(uold)*trans(id(v))*N())*idv(mesh_velocity)
                               +penalbc*trans(idv(*M_mesh_velocity))*id(v)/hFace() );

            }
#endif
        //--------------------------------------------------------------------------------------------------//

        //bool M_isStat=true;
        if (!M_isStat)
            {
                //termes temporels (ordre 1 pour le moment)
                form2( Xh, Xh, A ) +=
                    integrate( elements(mesh), _Q< 2*Order >(),
                               + trans(idt(u))*id(v)*M_bdf_fluid->polyDerivCoefficient(0));
                //element_fluid_type Buzz( Xh, "v" );
                //element_fluid_0_type buzz = U.template element<0>();
                //buzz.container() = M_bdf_fluid->polyDeriv().template element<0>().container();

                //auto buzz = M_bdf_fluid->polyDeriv().template element<0>();

                auto buzz = M_bdf_fluid->unknown(0).template element<0>();
                buzz.scale(M_bdf_fluid->polyDerivCoefficient(0));

                form1( Xh, F ) +=
                    integrate( elements(mesh), _Q< 2*Order >(),
                               + trans(idv(buzz))*id(v) );



            }

        //--------------------------------------------------------------------------------------------------//

#if STABILISATION

        //AUTO( beta_N , trans(N())*idv(beta) );
        AUTO( beta_abs , vf::sqrt( trans(idv(uold)-idv(*M_mesh_velocity))*(idv(uold)-idv(*M_mesh_velocity)) ));

        //terme de stabilisation
        form2( Xh, Xh, A ) +=
            integrate( internalfaces(mesh), _Q< 3*Order >(),
                       beta_abs*vf::pow(hFace(),2.0)*(1./vf::pow(cst(double(Order)),3.5))*trans(jumpt(gradt(u)))*jump(grad(v)) );

#endif

        //--------------------------------------------------------------------------------------------------//

        A->close();
        F->close();

        std::cout << "PtFixe : Assemblage finish\n";

        //--------------------------------------------------------------------------------------------------//

        if (!M_weakCL)
            {
                std::cout << "PtFixe : CL Fortes start\n";
                //form2( Xh, Xh, A ) +=
                //    on( markedfaces(mesh, mesh->markerName("Paroi")) , u, F, vec(cst(0.),cst(0.)) );
                //form2( Xh->template functionSpace<0>(), Xh->template functionSpace<0>(), A ) +=
                //on( markedfaces(mesh, mesh->markerName("Paroi")) , u, F, vec(cst(0.),cst(0.)),ON_ELIMINATION );
                form2( Xh, Xh, A ) +=
                    on( markedfaces(mesh, mesh->markerName("ParoiH")) , u, F, idv(*M_mesh_velocity) );//vec(cst(0.),cst(0.)) ); );
                form2( Xh, Xh, A ) +=
                    on( markedfaces(mesh, mesh->markerName("ParoiB")) , u, F, idv(*M_mesh_velocity) );//vec(cst(0.),cst(0.)) );//idv(mesh_velocity) );
                //form2( Xh, Xh, A ) +=
                    //on( markedfaces(mesh, mesh->markerName("Entree")) , u, F, Vit_Entree );
                    //on( markedfaces(mesh, mesh->markerName("Entree")) , p, F, vf::cst(10) );

                //element_fluid_velocity_type uYYY(Xh->template functionSpace<0>(),"xxx");
                //modifVec222(p,F, vf::cst(10),"Entree");
                //modifVec222(uYYY,F, vec(cst(0.),cst(0.)),"ParoiH");
                //modifVec222(uYYY,F, vec(cst(0.),cst(0.)),"ParoiB");


                std::cout << "PtFixe : CL forte finish\n";
            }

        //--------------------------------------------------------------------------------------------------//
        /*
        std::cout << "PtFixe : Solve start\n";
        vector_ptrtype X( M_backend->newVector( Xh ) );
        M_backend->solve( _matrix=A, _solution=X, _rhs=F,_maxit=1000 );
        U = *X;
        std::cout << "PtFixe : Solve finish\n";
    */

#endif
    }


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::updateStokes(const vector_ptrtype& X, sparse_matrix_ptrtype& A , vector_ptrtype& F)
    {

        using namespace Feel::vf;

        if (M_LookComment) std::cout <<"[NavierStokes] : updateStokes start \n";

        this->changeRepository();

        auto mesh = this->mesh();
        auto Xh = this->functionSpace();

        element_fluid_type U( Xh, "u" );U = *X;
        element_fluid_type V( Xh, "v" );
        element_fluid_0_type u = U.template element<0>();
        element_fluid_0_type v = V.template element<0>();
        element_fluid_1_type p = U.template element<1>();
        element_fluid_1_type q = V.template element<1>();


        //Tenseur des deformations
        auto deft = 0.5*(gradt(u)+trans(gradt(u)) );
        // auto def = 0.5*(grad(v)+trans(grad(v)) );

        //Matrice Identitee
        auto Id = oneX()*trans(oneX()) + oneY()*trans(oneY());
        // Tenseur des contraintes (trial)
        auto Sigmat = -idt(p)*Id + 2*M_mu*deft;

        // Tenseur des contraintes (test)
        //  auto Sigma = -id(q)*Id + 2*mu*def;

        //force volumique :
        //auto f = -rho*cst(9.81)*vec(cst(0.),cst(1.));
        auto f = vec(cst(0.),cst(0.));
        //auto f = FORCE_VOLUMIQUE ;

        form2( Xh,Xh, A, _init=true );

        // sigma : grad(v) sur Omega
        form2( Xh, Xh, A ) +=
            integrate( elements(mesh), _Q<2*(Order)>(),
                       + trace(Sigmat*trans(grad(v))) );

        // terme assurant l'incompressibilite
        form2( Xh, Xh, A ) +=
            integrate( elements(mesh),// _Q<2*(Order-1)>(),
                       + divt(u)*id(q) );

        // terme necessaire pour condition forte avec petsc
        form2( Xh, Xh, A ) +=
            integrate( elements( mesh ),// _Q<3*uOrder-1+3*(GeoOrder-1)>(),
                       +0*div(v)*idt(p)+0*divt(u)*id(q) +0*idt(p)*id(q) );

        // forces volumique
        form1( Xh, F, _init=true ) =
            integrate( elements(mesh), // _Q<Order>(),
                       trans(f)*id(v) );


        auto clDef = CL_DEF_NS(M_time);

        ForEachCL( clDef, cl::neumann_scal,
                   form1( Xh, F ) +=
                   integrate( markedfaces(mesh,mesh->markerName(PhysicalName)),// _Q<2*Order>(),
                              + trans( Expression*N())*id(v) ) );



        if (M_weakCL)
            {
                ForEachCL( clDef, cl::dirichlet_vec,
                           form2( Xh, Xh, A ) +=
                           integrate( markedfaces(mesh,mesh->markerName(PhysicalName)),/* _Q<2*Order>(),*/
                                      - trans(Sigmat*N())*id(v)
                                      + penalbc*trans(idt(u))*id(v)/hFace() ) );

                ForEachCL( clDef, cl::dirichlet_vec,
                           form1( Xh, F ) +=
                           integrate( markedfaces(mesh,mesh->markerName(PhysicalName)), /*_Q<2*Order>(),*/
                                      + penalbc*trans(Expression)*id(v)/hFace() ) );

                ForEachCL( clDef, cl::paroi_mobile,
                           form2( Xh, Xh, A ) +=
                           integrate( markedfaces(mesh,mesh->markerName(PhysicalName)),/* _Q<2*Order>(),*/
                                      - trans(Sigmat*N())*id(v)
                                      + penalbc*trans(idt(u))*id(v)/hFace() ) );

                ForEachCL( clDef, cl::paroi_mobile,
                           form1( Xh, F ) +=
                           integrate( markedfaces(mesh,mesh->markerName(PhysicalName)), /*_Q<2*Order>(),*/
                                      + penalbc*trans(idv(this->meshVelocity()))*id(v)/hFace() ) );

            }


        if (!M_isStat)
            {
                form2( Xh, Xh, A ) +=
                    integrate( elements(mesh), // _Q< 2*Order >(),
                               + rho*trans(idt(u))*id(v)*M_bdf_fluid->polyDerivCoefficient(0) );

                element_fluid_type Buzz( Xh, "v" );
                element_fluid_0_type buzz = Buzz.template element<0>();
                buzz.container() = M_bdf_fluid->polyDeriv().template element<0>().container();

                //auto buzz = M_bdf_fluid->unknown(0).template element<0>();
                //auto buzz = M_bdf_fluid->polyDeriv().template element<0>();

                form1( Xh, F ) +=
                    integrate( elements(mesh),// _Q< 2*Order >(),
                               - rho*trans(idv(buzz))*id(v) );
            }


        A->close();
        F->close();

        //std::cout << "Stokes : Assemblage finish\n";

        if (!M_weakCL)
            {
                ForEachCL( clDef, cl::paroi_mobile,
                           form2( Xh, Xh, A ) +=
                           on( markedfaces(mesh, mesh->markerName(PhysicalName)) , u, F, idv(this->meshVelocity()) ) );
                ForEachCL( clDef, cl::dirichlet_vec,
                           form2( Xh, Xh, A ) +=
                           on( markedfaces(mesh, mesh->markerName(PhysicalName)) , u, F, Expression ) );
            }

        if (M_LookComment) std::cout <<"[NavierStokes] : updateStokes finish \n";

    }

} // end namespace FSI
} // end namespace Feel


    template class Feel::FSI::NavierStokes<2,2,Feel::Simplex>;
//template class Feel::NavierStokes<2,4,Feel::Simplex>;
